A Quick Tour Through State Space reconstruction
Short Course ‘An Introduction to Empirical Dynamics Modelling, ICTP SAIRF’, Tutorial #1
This tutorial presents the general idea of atractor reconstruction with syntetic data, and there is a simulated and real data analysis proposed as exercises. For these activities, you will need the most recent version of R and the rEDM package installed in your work computer. To start, open R and load the package by typing the follwoing code in the R console:
library(rEDM)Before you proceed please read pass through the essential concepts by reading the introduction and the first section (“Empirical Dynamic Modelling”) of rEDM’s tutorial that you can find here or in your local R installation with the command
vignette("rEDM-tutorial", package="rEDM")A simple case: harmonic oscilator
This first part is to get the grips with R language and to figure out how to reconstruct a state space with only one time series and its lags. For the first example, we use a very well known system and its representation on state space, the bidimensional harmonic oscilator:
\(\frac{dx}{dt}=y\)
\(\frac{dy}{dt}=-x\)
which has the solution \(x(t)=\sin(t), \ y(t)=\cos(t)\) and so we can draw the state space of this system by laying down in two cartesian axis its two degrees of freedom, \(x(t)\) and \(y(t)\). This system has a circular atractor, that is, all the possible states lie over a circle in the state space:
## Generate the data at tau = 0.2 time intervals
time1<-seq(0, 500, by = 0.2)
x1<-sin(time1)
y1<-cos(time1)
## Plot the state space
plot(x1, y1, type = 'l', main = "State-space of harmonic oscilator",
cex=0.2, xlab="x(t)", ylab="y(t)")And here are the time series for the two state variables:
## Time series of the 1st 250 observations
ind <- 1:250
plot(time1[ind], x1[ind], type = 'l',
main = "Time series of the variables of the system",
col = "red", xlab = "Time", ylab = "Amplitude",
ylim=c(-1,1.25))
lines(time1[ind], y1[ind], col = "blue")
legend("topright", c("x","y"), col=c("red", "blue"),
lty=1, bty="n")A reconstruction of the state space
Now suppose that all you information you have about this system is a time series of one of the state variables, say \(x\). Because \(y(t)\) is missing we are not able to construct a state space as we saw above, but by using the result of Takens theorem we can construct a shadow manifold that recovers the information of the true state space. This reconstruction is not a copy of the original atractor, but it is a valid topological reconstruction of the original, a “shadow”.
To do this, we generate new variables by lagging the time series by multiples of a time interval \(\tau\). In this case we will use the interval between observations (\(\tau=0.2\) time units). The R commands below creates a data table object ( which is called a data frame in R) with the original time series and its lagged version. The handy function make_block by the rEDM authors does this job. Use this link to download the function to your working directory and then load it with the command:
source("../utilities/make_block.R")And now you can create a lagged variable dataframe in R by copying the following code and pasting in the R console:
## Data frame with observations with tau = 0.2 with one lag
df1 <- make_block(data.frame(x=x1), t=time1, max_lag = 2)A quick look at the first and last lines of the data frame can help to understand the lagged variabel structure:
## Take a look in the lagged data to understand how the function works
head(df1)## time x x_1
## 1 0.0 0.0000000 NA
## 2 0.2 0.1986693 0.0000000
## 3 0.4 0.3894183 0.1986693
## 4 0.6 0.5646425 0.3894183
## 5 0.8 0.7173561 0.5646425
## 6 1.0 0.8414710 0.7173561
tail(df1)## time x x_1
## 2496 499.0 0.49099533 0.65428133
## 2497 499.2 0.30813490 0.49099533
## 2498 499.4 0.11299011 0.30813490
## 2499 499.6 -0.08665925 0.11299011
## 2500 499.8 -0.28285377 -0.08665925
## 2501 500.0 -0.46777181 -0.28285377
By ploting the time series as a function of their lagged versions we get the shadow manifold, which is similar to the original one:
plot(x_1 ~ x, data=df1, type="l",
xlab = "X", ylab = "X(t + tau)" )There is a minimum number of dimensions to build a valid rceonstruction of the original atractor. In this simple case we obviously need two dimensions. Takens Theorem guarantees that the number of dimensions is up to 2D + 1, where D is the dimmensionality of the the attractor. Finding the optimal dimension for the shadow manifold is part of a vast topic called Embedology (Sauer et al. 1991). The package rEDM has some methods to find the optimal dimensions of the manifold. We will learn more in the following tutorials of this course.
Noisy data
How robust is the atractor reconstruction to random noise? There are two main types of noise or error in data, which we call process errors and observational errors. The first one is about random changes of the state variables, like environmental fluctuations that affect the dynamic. The second type of error appears when we are not able to measure the data with total precision, when we miss some counts on the counting of individuals of species, or do error in the sampling or measurements. Observational errors do not affect the dynamic, only the portrait we make of it.
In this section, we check how robust the atractor is to error by adding simulated noise to the original time series. Use the commands below to simulate that the time series has independent measurement errors with constant variance. To do this, we add to each observation a value drawn from a Gaussian distribution with zero mean and standard deviation which is 1/12 of the standard deviation of the time series itself.
## a new vector of observed values + Gaussian iid error
x2 <- x1 + rnorm(n = length(x1), mean = 0, sd = sd(x1)/12)
## Time series
plot(x1[ind] ~ time1[ind], type="l", xlab = "X", ylab = "Time",
col="grey", lwd = 5)
lines(x2[ind] ~ time1[ind], type="l", col="red", lwd=2)We thus proceed by creating a new data frame with the lagged variables and by ploting the new shadow manifold:
## Data frame with lagged variables
df2 <- make_block( data.frame(x = x2))
## Shadow manifold with embedding = 2
plot(x_1 ~ x, data=df2, type="l",
ylab = "X(t + tau)", xlab = "X" , col="grey",
lwd=0.5)In this case the shadow manifold still looks like the original shape of the attractor, despite the noise caused by measurement errors.
Increasing the noise
Measurement errors make atractors more complex by spreading points (and thus trajectories) over the state space. Hence, the trajectories may cross in the shadow manifold, which breaks down the one to one mapping to the original atractor. Of course this effect gets larger as we increase the measurement error. We now will investigate a bit more how this affects the atractor reconstruction. The commands below doubles the measurement error of the time series and plots the shadow manifold:
x3 <- x1 + rnorm(n = length(x1), mean = 0, sd = sd(x1)/6)
df3 <- make_block(data.frame(x=x3), max_lag = 5)
plot(x_1 ~ x, data=df3, type="l", col="grey", lwd=0.5,
ylab = "X(t + tau)", xlab = "X" )The atractor reconstruction looks worse now, but we can overcome this problem by adding more dimensions to the shadow manifold. By doing this we unfold the shadow manifold, that is unfolding crossing trajectories. We do this adding more axes, that correspond to the time series lagged by the further time steps (\(t + \tau, \ t + 2\tau, \ t + 3\tau\)…). Even the if the original manifold has a lower dimension its is valid to reconstruct the state space with more dimensions[^1]:
You can check with the interactive plot below how an additional dimension unfolds a bit more the reconstructed atractor. In the following days we will learn some methods to find the optimal dimensions of the manifold. By now the take-home message is: the more complex the atractor is, the more dimensions you will need to have a good reconstruction.
Problems
What is the effect of \(\tau\) ?
The commands below simulate a new time series of the bimensional harmonic oscilator without measurement error but with a much smaller interval between observations (\(\tau=0.01\) time units).
## Generate the data at 0.01 time intervals
time4<-seq(0, 500, by = 0.01)
x4<-sin(time4)
y4<-cos(time4)This change in \(\tau\) affects the shape of the reconstruction:
df4 <- make_block(data.frame(x=x4), t=time4)
plot(x_1 ~ x, data=df4, type="l",
xlab = "X", ylab = "X(t + tau)" )How do you explain that? Which practical issues does this fact bring and how to solve them?
Hints
1. The command above adds a line \(x(t) = x(t + tau)\) to the manifold, which can help you to understand the problem with small values of \(\tau\):
abline(0,1, col="red")2. Check this site: (http://www.physics.emory.edu/faculty/weeks//research/tseries4.html)
Atractor reconstruction
By using only the functions learned at this tutorial, try to reconstruct the state space of the time series] in this csv file, checking if a three-dimensional manifold improve the reconstruction. To load the data in R download the csv file and run the command:
prob2 <- read.csv("problem-time-series.csv")And to take look in the time series plot run:
plot( X ~ time, data = prob2, type="l")Hints
1. Keep in mind the insight from the previous problem: the value of \(\tau\) change the shape of the shadow manifold.
2. To create a simple 3D plot in R you can use the scatterplot3d package. The example below does this for one of the noisy timeseries we examined above:
## Loads the package, if not available
library(scatterplot3d)
scatterplot3d(x = df2$x, y = df2$x_1, z = df2$x_2,
type = "l", color="grey", lwd = 0.5,
xlab = "x(t)", ylab = "x(t+tau)", zlab = "x(t + 2 tau)")If the package is not available in your local R installation you can install it with:
install.packages("scatterplot3d")Learn more
- Sugihara, G., May, R., Ye, H., Hsieh, C.H., Deyle, E., Fogarty, M. and Munch, S., 2012. Detecting causality in complex ecosystems. science, 338(6106), pp.496-500.
- DeAngelis, D.L. and Yurek, S., 2015. Equation-free modeling unravels the behavior of complex ecological systems. Proceedings of the National Academy of Sciences, 112(13), pp.3856-3857.
- Taken’s Theorem, a video by Sugihara’s Lab.
- Takens, Floris. Detecting strange attractors in turbulence. Lecture notes in mathematics 898.1 (1981): 366-381.
- Deyle, E.R. and Sugihara, G., 2011. Generalized theorems for nonlinear state space reconstruction. PLoS One, 6(3), p.e18295.
- Sauer, T., Yorke, J.A. and Casdagli, M., 1991. Embedology. Journal of statistical Physics, 65(3), pp.579-616.
[^1] Sauer, Tim, James A. Yorke, and Martin Casdagli. “Embedology.” Journal of statistical Physics 65.3 (1991): 579-616.